/**
 * @author Ahammer
 * @date    2009 06
 * @update  2011 01 changed high[b]=typeGreyMax to high[b]=typeGreyMax+1;
 */

/*
 * This file is part of Iqm.
 * Copyright (c) 2010-2011 Helmut Ahammer
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package op;

import java.awt.image.renderable.ParameterBlock;
import javax.media.jai.Histogram;
import javax.media.jai.JAI;
import javax.media.jai.ParameterBlockJAI;
import javax.media.jai.PlanarImage;
import javax.media.jai.RenderedOp;
import main.Board;
import main.IqmTools;


/**
 * This is the main image processing class
 * There exist two approaches:
 * A user defined JAI operator is just called 
 * or:
 * The actual processing is implemented in this class 
 */
public class IqmHistoModifyOperator{ 


  public IqmHistoModifyOperator() {
	  //WARNING: Don't declare fields here
	  //Fields declared here aren't thread safe!   
  }
 
  private Histogram getHistogram (RenderedOp rendOp){
      Histogram histo = (Histogram) rendOp.getProperty("histogram");
      return histo;
  }
  private  RenderedOp getHistogramRendOp (PlanarImage pi){
	  int numBands = pi.getNumBands();
	  double typeGreyMax = IqmTools.getImgTypeGreyMax(pi);
      int[] bins    = new int[numBands];
      double[] low  = new double[numBands];
      double[] high = new double[numBands];
      for (int b = 0; b < numBands; b++){
    	  bins[b] = (int)typeGreyMax+1;
    	  low[b]  = 0.0d;
    	  high[b] = typeGreyMax+1;
      }
	  ParameterBlock pb = new ParameterBlock();	   
      pb.removeSources();
      pb.removeParameters();
      pb.addSource(pi);
      pb.add(null);
      pb.add(1);
      pb.add(1);
      pb.add(bins);
      pb.add(low);
      pb.add(high);
      RenderedOp rendOp = JAI.create("histogram", pb, null);
      return rendOp;
  }
  private int[][] getSubTotalValues(Histogram histo, double subLow, double subHigh ){
	  int numBands = histo.getNumBands();
	  int numBins  = histo.getNumBins(0);
      int[] totals = histo.getTotals();
      int[] subTotalLow  = new int[numBands];
      int[] subTotalHigh = new int[numBands];;
      int[][] subTotals  = new int[2][numBands]; // low and high
      for (int b = 0; b < numBands; b++){
    	  //calculate low subTotalLevel
    	  subTotalLow[b] = 0; //initial value
	      double subTotal = 0.0d;
	      while (subTotal/totals[b] < subLow){
	    	  subTotalLow[b] = subTotalLow[b] + 1; 
		      subTotal = (double)histo.getSubTotal(b, 0, (int)subTotalLow[b]);
	    	  if (subTotalLow[b] >= (numBins - 10)){
	    		  Board.appendTexln("IqmHistoModifyOperator: Alert: Low Offset value too high:");	
	    		  Board.appendTexln("IqmHistoModifyOperator: Low Offset value set to: " + subTotalLow[b]);
	    		  break;
	    	  }
		      //System.out.println("IqmHistoModifyOperator: subTotal: " +  subTotal);
	      }
//	      System.out.println("IqmHistoModifyOperator: subTotalLow[b]: " +  subTotalLow[b]);
//	      System.out.println("IqmHistoModifyOperator: subTotal/Totals[b]: " +  (subTotal/totals[b]));
	      //calculate high subTotalLevel
    	  subTotalHigh[b] = numBins-1; //initial value
	      subTotal = 0.0d;
	      while (subTotal/totals[b] < subHigh){
	    	  subTotalHigh[b] = subTotalHigh[b] - 1; 
		      subTotal = (double)histo.getSubTotal(b, (int)subTotalHigh[b], numBins-1 );
	    	  if (subTotalHigh[b] <= (subTotalLow[b] + 10)){
	    		  Board.appendTexln("IqmHistoModifyOperator: Alert: High Offset value too high:");	
	    		  Board.appendTexln("IqmHistoModifyOperator: High Offset value set to: " + subTotalHigh[b]);
	    		  break;
	    	  }
		      //System.out.println("IqmHistoModifyOperator: subTotal: " +  subTotal);
	      }
//	      System.out.println("IqmHistoModifyOperator: subTotalHigh[b]: " +  subTotalHigh[b]);
//	      System.out.println("IqmHistoModifyOperator: subTotal/Totals[b]: " +  (subTotal/totals[b]));     
      }
      subTotals[0] = subTotalLow;
      subTotals[1] = subTotalHigh;     
	  return subTotals;
  }

  /**
   * 
   * @param ParametrBlockJAI pb
   */
  public Object run(ParameterBlockJAI pbJAI){
	  Object ob = null;
	  //ob = JAI.create("IqmHistoModify".toLowerCase(), pb, null);  
	  
	  PlanarImage pi  = (PlanarImage) pbJAI.getSource(0);
	  //IqmTools.plotHistogram(pi);
	  String imgName = (String) pi.getProperty("image_name");
	  double typeGreyMax = IqmTools.getImgTypeGreyMax(pi);
	  int numBands = pi.getNumBands();
	  int width    = pi.getWidth();
	  int height   = pi.getHeight();
	  int method   = pbJAI.getIntParameter("Method");
	  double offSetLow  = (double)pbJAI.getIntParameter("OffSetLow")/1000.0d; //in prozent
	  double offSetHigh = (double)pbJAI.getIntParameter("OffSetHigh")/1000.0d;  
	  boolean isBandCombineSet = false;
	  if (pbJAI.getIntParameter("BandCombine") == 0) isBandCombineSet = false;
	  if (pbJAI.getIntParameter("BandCombine") == 1) isBandCombineSet = true;


	  if ((offSetLow > 0)||(offSetHigh > 0)){
		  System.out.println("IqmHistoModifyOperator: offSetLow: " +  offSetLow);
		  System.out.println("IqmHistoModifyOperator: offSetHigh: " +  offSetHigh);
		  //IqmTools.plotHistogram(pi);

		  RenderedOp rendOp = getHistogramRendOp(pi);
		  Histogram histo = getHistogram(rendOp);

	      int[][] subTotals = getSubTotalValues(histo, offSetLow, offSetHigh);
	      
		  double[] lowClamp  = new double[numBands];  //= {offSetLowValue,  offSetLowValue,  offSetLowValue};
		  double[] highClamp = new double[numBands]; //= {offSetHighValue, offSetHighValue, offSetHighValue};
		  for (int b = 0; b < numBands; b++){
			  lowClamp[b]  = subTotals[0][b];
			  highClamp[b] = subTotals[1][b];
		      System.out.println("IqmHistoModifyOperator: Band, lowClamp[b], highClamp[b] "+b+" " +  lowClamp[b]+"  "+highClamp[b]);
		  }
		  if (isBandCombineSet){
			  double lowClampMean  = 0d;
			  double highClampMean = 0d;
			  for (int b= 0; b < numBands; b++){
				  lowClampMean  = lowClampMean  + lowClamp[b];
				  highClampMean = highClampMean + highClamp[b];
			  }
			  lowClampMean  = lowClampMean/numBands;  //Mean
			  highClampMean = highClampMean/numBands;
			  for (int b= 0; b < numBands; b++){
				  lowClamp[b]  = lowClampMean;
				  highClamp[b] = highClampMean;
			      System.out.println("IqmHistoModifyOperator: Band, lowClamp[b], highClamp[b] "+b+" " +  lowClamp[b]+"  "+highClamp[b]);
			  }
		  }
		  //IqmTools.plotHistogram(pi);
		  ParameterBlock pb = new ParameterBlock();	   
	      pb.removeParameters();
	  	  pb.removeSources(); 
	      pb.addSource(pi);  
	      pb.add(lowClamp);
	      pb.add(highClamp);
	      pi = JAI.create("clamp", pb, null);
		  //IqmTools.plotHistogram(pi);
	  }
	
	  if (method == 0){ //rescale
	      //get Extrema
	      ParameterBlock pb = new ParameterBlock();
	      pb.addSource(pi);
	      pb.add(null);
	      pb.add(1);
	      pb.add(1);
	      RenderedOp rendOp = JAI.create("extrema", pb, null);
	      double[] greyMinVec = (double[])rendOp.getProperty("minimum");  
	      double[] greyMaxVec = (double[])rendOp.getProperty("maximum");  
	      double greyMin = greyMinVec[0];  
	      double greyMax = greyMaxVec[0];  
	      for(int b = 1; b < greyMinVec.length; b++)  
	       	      {  
	       	      if (greyMinVec[b] < greyMin) greyMin = greyMinVec[b];  
	       	      if (greyMaxVec[b] > greyMax) greyMax = greyMaxVec[b];  
	      }
	      double[] scale  = new double[1];
	      double[] offset = new double[1];
	      scale[0]  = typeGreyMax/(greyMax-greyMin);  
	      offset[0] = typeGreyMax*greyMin/(greyMin-greyMax);  

	      pb.removeParameters();
	  	  pb.removeSources(); 
	      pb.addSource(pi);  
	      pb.add(scale);  
	      pb.add(offset);  
	      ob = (PlanarImage)JAI.create("rescale", pb, null);
	      //IqmTools.plotHistogram((PlanarImage) ob);
	  }
	  if (method == 1){ //equalize   //Flat uniform  histogram 
		  float[][] CDFeq = new float[numBands][];
		  for (int b = 0; b < numBands; b++){
			  CDFeq[b] = new float[(int)typeGreyMax+1];
			  for (int i = 0; i <= typeGreyMax; i++){
				  CDFeq[b][i] = (float)(i+1)/((float)(typeGreyMax+1));
			  }
		  }
		
		  RenderedOp rendOp = getHistogramRendOp(pi);
	      ob = JAI.create("matchcdf", rendOp, CDFeq);
	      //IqmTools.plotHistogram((PlanarImage) ob);
	      
		  //old method
//	      int sum = 0;
//	      byte[] cumulative = new byte[(int)typeGreyMax+1];
//	      int array[] = getHistogram(pi, typeGreyMax);
//	      float scale = (float) (typeGreyMax / (float) (width * height));           
//	      for ( int i = 0; i <= typeGreyMax; i++ ) {
//	           sum += array[i];
//	           cumulative[i] = (byte)((sum * scale) + .5F);
//	      }
//	      LookupTableJAI lookup = new LookupTableJAI(cumulative);
//	      ParameterBlock pb = new ParameterBlock();
//	      pb.addSource(pi);
//	      pb.add(lookup);
//	      ob = JAI.create("lookup", pb, null);
//	      IqmTools.plotHistogram((PlanarImage) ob);
	  }
	  
	  if (method == 2){ //normalize
		  double[] mean  = new double[] { typeGreyMax/2d, typeGreyMax/2d, typeGreyMax/2d };
	      double[] stDev = new double[] { (typeGreyMax+1d)/4d, (typeGreyMax+1d)/4d, (typeGreyMax+1d)/4d  };  //?? 34 bei Byte  bei short???
	      float[][] CDFnorm = new float[numBands][];
	      for (int b = 0; b < numBands; b++){
		      CDFnorm[b] = new float[(int)typeGreyMax+1];
		      double mu = mean[b];
		      double twoSigmaSquared = 2.0*stDev[b]*stDev[b];
		      CDFnorm[b][0] = (float)Math.exp(-mu*mu/twoSigmaSquared);
		      for ( int i = 1; i <= typeGreyMax; i++ ) {
		            double deviation = i - mu;
		            CDFnorm[b][i] = CDFnorm[b][i-1] + (float)Math.exp(-deviation*deviation/twoSigmaSquared);
		      }
	      }
	      for (int b = 0; b < numBands; b++){
		      double CDFnormLast = CDFnorm[b][(int)typeGreyMax];
		      for ( int i = 0; i <= typeGreyMax; i++ ) {
		            CDFnorm[b][i] /= CDFnormLast;
		      }
	      }
    
		  RenderedOp rendOp = getHistogramRendOp(pi);
	      ob = JAI.create("matchcdf", rendOp, CDFnorm);
	      //IqmTools.plotHistogram((PlanarImage) ob);
	  }
	  if (method == 3){ //piecewise
		  float[][][] bp = new float[1][2][];
	      bp[0][0] = new float[] { 0.0F,  32.0f/255.0f*((float)typeGreyMax),  64.0f/255.0f*((float)typeGreyMax), (float)typeGreyMax };
	      bp[0][1] = new float[] { 0.0F, 128.0f/255.0f*((float)typeGreyMax), 112.0f/255.0f*((float)typeGreyMax), (float)typeGreyMax };
	      //bp[0][1] = new float[] { 0.0F, ((float)typeGreyMax+1f)/2f, 112.0f, (float)typeGreyMax };
	      //bp[0][0] = new float[] { 0.0F, 32.0F, 64.0F, 255.0F };
	      //bp[0][1] = new float[] { 0.0F, 128.0F, 112.0F, 255.0F };
	      ob = JAI.create("piecewise", pi, bp);
	      //IqmTools.plotHistogram((PlanarImage) ob);
	  }
	  
	  //((PlanarImage)ob).setProperty("image_name", imgName);
	  return ob;
  }
}
